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ABSTRACT 

We develop and calibrate a realistic model flame for hydrodynamical simulations of deflagrations in 
white dwarf (Type la) supernovae. Our flame model builds on the advection-diffusion-reaction model 
of Khokhlov and includes electron screening and Coulomb corrections to the equation of state in a 
self-consistent way. We calibrate this model flame — its energetics and timescales for energy release 
and neutronization — with self-heating reaction network calculations that include both these Coulomb 
effects and up-to-date weak interactions. The burned material evolves post-flame due to both weak 
interactions and hydrodynamic changes in density and temperature. We develop a scheme to follow 
the evolution, including neutronization, of the NSE state subsequent to the passage of the flame front. 
As a result, our model flame is suitable for deflagration simulations over a wide range of initial central 
densities and can track the temperature and electron fraction of the burned material through the 
explosion and into the expansion of the ejecta. 

Subject headings: hydrodynamics — nuclear reactions, nucleosynthesis, abundances — supernovae: 
general — white dwarfs 



1. INTRODUCTION 

Type la supernovae (SNe) are bright explosions charac- 
terized by strong silicon P Cygni features near maximum 
light and a lack of hydrogen spectral features. The cur- 
rently favored interpretation is the disruption of a near- 
Chandrasekhar-mass C/0 whit e dwarf by a thermonu- 
clear runaway (for a review, see iHillebrandt fc Niemeverl 
120001 and references therein). These events are fasci- 
nating in and of themselves and are important both for 
their contribution to the cosmic abundance of iron-peak 
elements and for their role as standard candles. 

Models of Type la SNe necessarily involve a mechanism 
for incinerating the star by a thermonuclear runaway, and 
the nature of this mechanism is the subject of contem- 
porary research. In the explosion, a thermonuclear flame 
propagates through the C/0 fuel o f the white dwarf as 
eithe r a subsonic deflagrat i on front (iNomoto et al]ll976l 
IT981 iReinecke et al.ll200l iGamezq et al.ll2003D or a su- 



personic detonation wave (|Arnettlll969t feoisseau et all 
Hogs') and releases sufRcient energy to unbind the star. 
However, models involving either a pure deflagration 
or a pure detonation have traditionally been unable 
to provide an explanation for both the observed ex- 
pansion velocities and the spectra produced by ejecta 
that are rich in intermediate-mass and iron-peak el- 
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ements ("Tnu-an & Canicron| Il971t iWooslev & Weaveil 
[1986: Petschek . Recent work, however, suggests 

that a fast deflagration alone may provide sufficient en- 
ergy to unbind the star ijHillebrandt fc Roepkell20d5|) . 

There has been considerable progress recently in hy- 
drodyn amical si mulations of deflagration s of C/0 white 
dwarfs ()R6pke fc H illebrandt 2005: Gamezo et al.ll2005(l 
that model the entire star. This is a complicated en- 
deavor, predominantly due to the vast range of length 
scales: the laminar flame width is ~ 10~-'-10cm, some 8 
to 12 orders of magnitude smaller than the stellar radius 
Irimmes fc Wooslcy 1993). Because the computational 
requirements for simulations with these disparate scales 
demand resources well beyond current capabilities, mul- 
tidimensional Type la models must make use of an ap- 
propriate sub-grid-scale model for the evolution of the 
thermonuclear burning front. Moreover, large-scale sim- 
ulations are very demanding of computational resources, 
and it is not feasible at present to include enough nuclides 
to allow for directly computing the reaction kinetics. A 
realistic model must accurately describe the nuclear en- 
ergy that is released, the timescale on which it is released, 
and the compositional changes that occur in the flame. 
In addition, the burned material continues to evolve after 
the passage of the flame due to both weak interactions 
and hydrodynamic evolution, and realistic simulations 
must describe this "post-flamc" evolution. 

In this paper we present a study of the nuclear burn- 
ing that occurs during C/0 deflagrations, with the goal 
of producing a realistic flame model for simulations of 
Type la supernovae. Building on the advection-diffusion- 
reaction (ADR) flame model of|^hokhlov (1995), we fur- 
ther develop a three-stage flame model. The work de- 
scribed here improves upon earlier calculations in sev- 
eral ways. First, screening and the Coulomb interaction 
are included self-consistently in the forward and inverse 
rates, so that detailed balance is preserved. Second, we 
account for the neutronization and neutrino loss rates 
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via electron captures at high densities. Our treatment is 
sufficiently general that it can also describe the expan- 
sion following the explosion. With minimal tuning, the 
method should capture the bulk energetics and associ- 
ated "freezing" of the abundance pattern as the matter 
expands. 

The present paper is concerned with describing with 
the utmost care the nuclear processes that occur in the 
real flame and proposing a method that can capture the 
necessary features: energy release, its timescales and 
neutronization. The hydrodynamic behavior in the 3- 
dimensional flow (importance of flame front curvature, 
acoustic behavior, flame front stability, effects of finite 
resolution, etc.) is under active investigation. Due to 
the complexity of those issues, the fact that the work 
is still in progress, and the clean separation of much of 
them from the nuclear physics addressed here, we have 
chosen to publish this work in advance of a forthcoming 
separate hydrodynamical study. It is important to note 
that such a hydrodynamical study must be performed in 
order to understand how this particular flame-capturing 
technique, with realistic energetics, behaves in simula- 
tions of a deflagration in a WD. 

Our three-stage flame model describes the propagation 
of a nuclear flame through a imiform '^C/'^O mixture. 
In this work we consider the case of a mixture with 1:1 
mass fractions, but t he method can b e applied to any 
ratio. As described in iKhokhlovl l)2000|) . the burning oc- 
curs in roughly three stages with well-defined timescales. 
First is the '^C-l-'^C fusion, leading to a mixture consist- 
ing of the unburned '^O together with ^"Ne, ^"^Mg, and 
a-particles. Second, the resulting mix burns on a longer 
timescale leading to the formation of predominantly Si- 
group (intermediate mass) a-elements and a-particles, a 
state that is commonly referred to as nuclear statistical 
quasi-equilibrium (NSQE). Third, these elements burn 
on a still longer timescale to nuclear statistical equilib- 
rium (NSE), which at the densities of interest in the de- 
flagration phase of Type la SNe, consists primarily of 
Fe-peak nuclei, a-particles, and protons. 

After the flame has burned the C/0 fuel, two factors 
continue to influence the energetics: (1) decreases (or 
possibly increases) in temperature and density resulting 
from the large-scale motion of the star, and (2) neutron- 
ization of the hot nuclear ash. Neutronization can heat or 
cool the ash, depending on whether the nuclear binding 
energy released by neutronization exceeds or falls short 
of the energy lost by neutrinos, and the competition be- 
tween these two processes is a sensitive function of den- 
sity. We note that because the light curves of SNe la 
are strongly sensitive to the mass of ^^Ni produced, the 
degree to which neutronization reduces the mean ^^Ni 
concentration over the regions of the core (extending out 
to « 0.8-1 M0) in which conditions of nuclear statisti- 
cal equilibrium are achieved is clearly a critical issue. It 
is therefore important to provide an accurate measure of 
the influence of electron captures on the energetics. Only 
in this way can we ensure that the peak temperatures 
and temperature histories of fluid elements are accurate. 
Accordingly, our flame model follows this neutronization. 

In this paper we describe nuclear flame propagation 
(§ 12)). We describe fully-resolved simulations of ther- 
monuclear flames propagating through a slab of stellar 
material and compare these to similar "one-zone" sim- 



ulations of nuclear burning. The "one-zone" or "self- 
heating" simulations allow for a longer-time evolution 
than the fully-resolved flames and provide the energet- 
ics and timescales for the model flame. Comparison to 
the fully-resolved flame calculations is an important con- 
sistency test. In this section we also describe the method 
by which we determine the density contrast across the 
flame as is required as input by the flame model and 
present contrasts at selected densities measured from the 
fully-resolved flame simulations. In (§ we describe 
our model flame. We outline the basics of our ADR 
flame capturing scheme and describe the stages of nu- 
clear burning and the energy release occurring in each. 
In (§^, we describe the calculations of energetics and 
timescales that serve as input to the flame model. We 
present the results of self-heating simulations of nuclear 
burning with a detailed nuclear reaction network and 
contemporary reaction rates. We report on the effect of 
including Coulomb corrections and screening and present 
timescales for burning to NSQE and NSE. In (§EIl, we 
describe the change in energy resulting from the evolu- 
tion of the burned material in NSE. We present a study 
of the NSE state that quantifies the change in binding 
energy with temperature and density evolution and the 
effects of weak interaction (neutronization) . We also de- 
scribe our method of representing the NSE state with 
a reduced set of nuclei. In (§ ^ we present results of 
model flame simulations at several densities, and in (§[7J) 
we draw conclusions from this effort. Finally, we present 
an appendix with details of the plasma coulomb correc- 
tions applied to both the network and NSE calculations. 

2. NUCLEAR FLAME PROPAGATION 

A thermonuclear flame in the interior of a star propa- 
gates subsonically. This slow propagation makes follow- 
ing the flame with an explicit hydrodynamics method 
such as PPM difficult because of the high spatial res- 
olution needed to resolve the flame. The time step 
in an explicit hydrodynamics method is limited by the 
Courant condition; namely, the code must resolve the 
sound-crossing time of the smallest zone. Because of this 
constraint, the time steps for a high- resolution simulation 
are very small and a large number of time steps is there- 
fore required to propagate the subsonic flame, even across 
just a few simulation zones. This problem is particularly 
acute at lower densities for which the timescale for burn- 
ing is relatively long. A fully-resolved flame simulation 
could literally take billions of time steps, which is imprac- 
tical and would most likely produce inaccurate results. 
We implement the fla me model in the explicit hydro- 
dyna mics code FLASH ijFrvxell et alJ 1200(1: ICalder et alJ 
|2002|) . but rely for input to the model flame on self- 
heating network calculations described in this section. 

2.1. Nuclear Flame Structure and Self-Heating 
Calculations 

Because of the difficulty encountered with propagat- 
ing flames with flash, particularly at low densities, we 
made use of self-heating network calculations (see section 
01 for details). These are (energetically) closed-box cal- 
culations of nuclear burning at either constant pressure 
(isobaric) or constant density (isochoric). We use both 
cases in this paper, each for a different purpose. Isobaric 
behavior should be quite close to what occurs in a fluid 
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Fig. 1. — Top panel: Results from a fully resolved DNS of a prop- 
agating flame at lO' g cm"' performed with FLASH and the aproxlO 
network. Shown are abundances for p, "^He, '^C, '*0, ^"Ne, ""^Mg, 
^*Si, ^''Fe, and ^'Ni plotted against the distance from the flame po- 
sition so that the flame is propagating to the left. The small "tails" 
at the right are artifacts of the "match head" used to ignite the 
flame. The simulation was run on an adaptive mesh with an ef- 
fective resolution of 7.81 X lO"^" 
n9 



"cm. Bottom panel: Results from an 
isochoric simulation at 10'' g cm"' performed with the self-heating 



19 nuclide network aproxl9. Shown are the same abundances as 
the top panel plotted against the time and the calculated distance. 
The distance scaling along the top axis was obtained by multiply- 
ing the time by the flame speed with respect to the (expanded) 
ash. 



element as the subsonic flame sweeps through it, so that, 
in a sense, such a sclf-hcating network calculation may be 
thought of as a "one-zone" calculation. Isochoric calcula- 
tions arc thermodynamically simpler and for this reason 
are used in section^to gauge the effects of screening and 
to compare different reaction networks. As will be shown 
below, the results of isobaric and isochoric calculations 
are fairly similar for high density (p > 10^ g cm"-'). 

We refer to an actual propagating flame simulated with 
FLASH using thermal diffusion as a direct numerical sim- 
ulation (DNS). An important test is the consistency be- 
tween actual flames propagated with flash and the re- 
sults of the self-heating calculations. We performed sim- 
ulations of propagating flames with flash at several rel- 
atively high densities for which a flame calculation was 
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an equal-mass mix 
"' Shown are the 



of C and O from an initial density of lO' g cm 
results from a isochoric {solid lines) and an isobaric (/" = 5 X 10^* 
dyne cm"^; dashed lines) self- heating calculation. Our 200 nuclide 
network was used in both calculations. The approximate distance 
scale is as in Fig.lTI 



feasible. The most notable difference between the DNS 
and the self-heating calculations is that it is necessary 
to use a fairly small reaction network (19 nuclides) in 
the DNS, whereas the self- heating calculations are able 
to use a much larger network (200 nuclides, see section 
EI). The DNS employed adapt ive mesh refinement with 
an effective resolution of 7.81 x 10"*'cm. The minimum 
refinement of the mesh was kept at 1.56 x 10"^ cm, and 
one additional level of AMR allowed the extra resolution 
of the flame front. 

A comparison of a DNS and isochoric self-heating is 
shown in Fig. ^ Both of these calculat i ons u s e the 19 
nucli de network aproxl9 l) Weaver et al.l 119781 iTimmed 
|1999() to burn fuel which is initially an equal mass mix- 
ture of '^C and '^O at a density of 10^ g cm"-'. The top 
panel shows how abundances vary with distance behind 
the flame front in the DNS (so that unburned material 
is on the left). The flame was ignited at the edge of a 
0.08 cm domain by raising the temperature of a small 
region of the domain to the unphysically high value of 
lO'^K. Heat then diffuses from this preheated region and 
ignites the flame, which then propagates across the grid. 
The rightmost parts of the mass fraction curves, at the 
largest distance from the flame, have small "tails" that re- 
sult from the unphysical initial conditions of the very hot 
"match head." Shown are mass fractions for p, ^He, '^C, 

16q 20jyg^ 24]y^g^ ISgj^ 54pg^ ^^^^ 56p^j^ rp^^ ^^^^^ ^^^^^ 

calculated from this simulat ion, 34.2 km s"^, agreed well 
with the resuh of Timnics fc WooslevI l|1992|) . 

The abundance evolution with time from an isochoric 
sclf-heating calculation using the same network is shown 
in the bottom panel of Fig. ^ We have calculated an 
approximate distance from the flame front (shown on 
the upper scale) by multiplying the time by the flame 
speed with respect to the ash. From conservation of mass 
this is related to the laminar flame speed in the fuel by 

5ash = 5lam(l+A)/(l - A), wllCrC A = (pfuel " Pash)/(Pfuel + Pash) 
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is the Atwood number. For this calculation we used val- 
ues of Ap2oo in Table ^ fse e description in S 12.211 an d 
the laminar flame speed from lTimmes &: WooslevI (|1992|) . 
5iam = 5.3 X 10'' cms-i for p = 10^ g cm'^ and XC^C) = 0.5 
(see eq. [H]). We find good agreement in terms of abun- 
dance structure between the isochoric self-heating calcu- 
lation and DNS. The stretching on the left in the bot- 
tom panel is due to the choice of zero point. With an 
initial temperature T = 1.7 x 10^ K, there is a long smol- 
dering phase in the self-heating calculation before rapid 
consumption of '^C occurs. The same "pre"-heating is 
accomplished in the DNS by thermal dilfusion, giving a 
heating time and length scale similar to the '^C consump- 
tion scales. To be able to broaden the region where C/0 
decreases sharply, we chose the time zero to be a time 
where C/0 starts to decrease. Since the C consumption 
timescale is short compared to the later flame stages (e.g. 
reaching NSE), this choice does not affect the total flame 
thickness. Also the depletions of '^C and '^O are shifted 
slightly further in distance in the self-heating plot by our 
use of a constant distance scaling because this is in reality 
where the expansion takes place. 

The next step is to compare these results with those 
of the much larger nuclear network that will be used for 
the rest of our self-heating calculations. Fig. |21 shows 
the abundance evolution for p = 10^ g cm~^ using the 
200 nuclide network described in § 0] The two networks 
give very similar evolutions, with the most notable dif- 
ference being in the final abundances. The 200 nuclide 
network has many more nuclides available near the Fe 
peak, such that less material is concentrated in a single 
nuclide (^"^Fe) than with the small network. The 5 most 
abundant nuclides in the NSE state are ^'*Fe, "^He, ^^Co, 
s^Ni, 5^Co with abundances of 16%, 13%, 8.4%, 7.3% and 
6.2% respectively, and there are another 14 nuclides with 
abundances between 1 and 5%. 

The evolution of a fluid element as the subsonic flame 
passes is in fact nearly isobaric. Due to the high level 
of degeneracy a.t p= lO' g cm~^, however, there is little 
expansion and so the isochoric and isobaric self-heating 
calculations give quite similar results. In Fig. |21 we also 
show the result from an isobaric calculation. Here the 
pressure (P = 5 x 10^^ dynescm"^) is chosen to give an 
initial density of lO' g cm"^. We can see that the iso- 
baric calculation takes slightly longer to reach NSQE, a 
result of a slightly decreased temperature due to expan- 
sion. The density is decreased to 8.3 x lO' g cm"-' at this 
phase. However, the NSE timescale and the final NSE 
abundances don't differ much from the isochoric calcula- 
tion. 

Before discussing the general features of this flame 
structure, it is useful to see how it changes when burn- 
ing at lower densities. Figs. 13 and 01 show isobaric self- 
heating calculations for initial densities of pi =10** and 10^ 
g cm"-'. Light particles, particularly a-particles, become 
less dominant in the NSE flnal state. It is already evi- 
dent that at these low densities the flame becomes very 
extended physically, such that a proflle like that shown 
in Fig. ^ could never be realized in an actual star. The 
dynamical response time of the star is shorter than these 
timcscales (see Fig. Ill|l , causing the reactions to freeze 
out due to expansion, and the scale height at this density 
is smaller than the large scales of the structure shown. 
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Fig. 3. — Similar to Fig. I^but for isobaric only with an initial 
density of 10* g cm"^. 
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Fig. 4. — Similar to Fig. |^ but for isobaric only with an initial 
density of lO' g cm"^. 



At all densities, though the overall progression 
timescale differs significantly, the three burning stages 
are evident. Beginning with equal parts by mass of '^C 
and '^O, initially the burning of the '~C leads to a mix- 
ture of '^O, ^"Nc, ^"^Mg, and a-particles. This mixture 
then burns to intermediate mass (^^Si to '^''Ca) nuclei and 
a-particles (NSQE) in which most of the energy release 
has taken place. Finally, the mixture reaches full NSE 
consisting of a temperature and density dependent mix- 
ture of iron-peak nuclei, a-particles, and protons. The 
transition points at p = 10^ g cm"^ (Fig- 01) ^.rc at 0.2, 
0.4, and 30 cm on the approximate distance scale. While 
the '^C burning and the transition to NSQE are well sep- 
arated at the low density, at the higher densities these 
stages appear more as an offset between consumption of 
i^C and i<'0. 

It is important to note again that the NSE state is not 
static. The forward and inverse rates are large but bal- 
anced, which enables the NSE state to adjust as the fluid 
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element evolves under hydrodynamical motions (see § |S1 
for further details). The flame model must reproduce this 
behavior with good approximations to the flame propa- 
gation speed, energy release timcscales, and total energy 
release. 

2.2. Inferring Density Contrast From Isochoric 
Calculations 

In the flamclet regime, the flame is accelerated by 
the Rayleigh- Taylor instability, which increases the sur- 
face area of the fla me and its effective propa gation 
speed (|Khokhlovlll99,^lNiemever fc Kersteinill997^ . The 
strength of this effect is related to the density contrast, 
characterized via the Atwood number. A, due to the en- 
ergy release in nuclear burning. Energy is released and 
does work to expand the matter and raises the tempera- 
ture until NSE is reached (for densities ^ 10^ g cm"-') or 
the nuclear energy is exhausted upon fusion to ^^Ni. 

The energy release found from an isochoric sclf-heating 
calculation can be used to estimate the expansion ex- 
pected for the actual (approximately isobaric) burn. We 
denote this approximate Atwood number by A, and cal- 
culate it by using a single-stage ADR flame (see Section 
that has a specified energy release, AQ, and ash 
composition taken from an isochoric self-heating calcu- 
lation (these are given in Table OJ. The ash composi- 
tion defines both the heat capacity and the relationship 
between the ion pressure contribution and the tempera- 
ture. Tabic ^ compares Atwood numbers calculated us- 
ing high-resolution DNS of flames, Adns, as was shown in 
Fig. n and A for isochoric self-heating calculations com- 
puted with two reaction networks, aproxlO and P200. By 
comparing these three cases we can see what differences 
are due to the nuclear network alone, and what are due 
to the isochoric approximation. Self-heating calculations 
with the aproxlQ nuclear network, the same as that used 
in the DNS, give a larger total energy release at the high 
densities relative to those found using our 200 nuclide 
network (see column 5), which is reflected in A (column 
3 and 4). 

The approximate Atwood numbers. A, are lower than 
the actual Atwood numbers. The ash temperatures, T^sh, 
found for the P200 case in this approximate scheme are 
shown in column 6, and are lower than the final temper- 
atures (shown in Table of the self heating calculations 
from which the AQ was drawn, reflecting the fact that 
some of the energy was used to do work for the expan- 
sion. This lower temperature means that the ash is now 
out of NSE, and if allowed to relax it would release nu- 
clear binding energy, expand, and raise its temperature 
to reach the A found in the DNS. 

We find that, at these densities (p > 10^ g cm"-'), iso- 
choric self-heating calculations provide reasonable esti- 
mates for the expansion factor and very gross estimates 
for the energy release. From this we conclude that A 
from the isochoric self-heating calculation can be used as 
a static approximation where needed. We should empha- 
size, however, that in the model flame described below, 
the actual total energy release of the flame is determined 
by the NSE in post-flame material, which is calculated 
dynamically as described in section [S] and therefore will 
reproduce that of an isobaric flame. 



3. IMPLEMENTATION OF MODEL FLAME ENERGETICS 

The key ingredient in numerical simulations of the de- 
flagration phase of SNc la is the nuclear flame model. In 
this section we describe our flame model, the implemen- 
tation of the energetics, and the implementation of the 
"post flame" NSE evolution. 

3.1. Basics of an ADR Flame model 

We us e a flame-ca pturing scheme based on the ADR 
m ethod llKhokhlovll 995. 200^ that is described in detail 
in lVladimirova et al.l H2006fl . The flame front is localized 
via the value of a reaction progress variable, </>, such that 
in the reactant = 0, in the product 0=1, and <f> varies 
monotonically across the flame front (see Fig. EJ. The 
width of the ADR flame is often much larger than the 
thickness of the physical flame it represents, as it is pre- 
scribed to be several computational zones thick, thereby 
becoming even kilometers in width in full-star supernova 
simulations. The value of (j) is defined everywhere on 
the grid, and can be associated with the mass fraction 
of burned material in a cell, provided that the reactant 
composition is initially uniform. 

The reaction progress variable is evolved via an 
advcction-diffusion- reaction equation, 

^ + vV0 = KV20+i7?((/)) , (1) 
at T 

with artificial reaction and diffusion coefficients 

where 0o marks the value of at which the reaction 
begins. 

When V = 0, the solution of equation is a traveling 
wave with speed So with a specific functional form, 4'{x— 
Sot). Traveling speed So and front thickness /q depend on 
the diffusivity, the reaction time, and the amplitude of 
the reaction rate, oc ^JRok/t and Iq oc kt/Ro- For 
the reaction rate lO, the expressions for So a.nd (l)(x-Scit) 
can be found analytically. As shown in lVladimirova et alJ 
l)2006j) . if the amplitude of the reaction rate Ro satisfies 
the relation {(f)Q-Ro)e^^'^'' +Ro = 0, the front propagates 
with the speed So = ^JkJt . In our implementation, 0o = 
0.3 and Rq = 0.3128, so the front speed and thickness obey 
the relations = \/k/t and to = ^-Jkt. 
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Fig. 5. — flame front structure for a one-stage ADR flame 
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TABLE 1 

Atwood numbers 



density 


j4dns 


Tash DNS 


■^aproxl9 


Ap200 




Tash 


(lO' g cm-3) 




10' g 






(10" erg g-l) 


(10' K) 


6 


0.065 


10.5 


0.056 


0.051 


0.21 


6.7 


2 


0.093 


9.34 


0.082 


0.077 


0.15 


6.8 


1 


0.116 


8.51 


0.101 


0.097 


0.08 


6.6 


0.5 


0.150 


7.80 


0.127 


0.125 


0.05 


6.3 


0.1 






0.227 


0.229 


-0.06 


5.3 



When V y 0, the front propagates with the speed 
with respect to the background velocity, provided the 
velocity variations on the scale of the front thickness are 
negligible in comparison with Sq. 

The ADR flame model uses the evolution of cj) accord- 
ing to equation ^ to propagate a front. The front speed 
S and the front thickness / are input parameters; they 
determine the coefficients k and r in equation (Q, 

K = fl5//4, T = //(4a5). (3) 
Here, a is an adjustment coefficient that accounts for 
the front speed dependence on thermal expansion (since 
thermal expansion results in velocity variations across 
the front and alters the front speed). The coefficient 
a depends on the density ratio a = Pu/ Pb and, for the 
reaction rate (El, was estimated as a = 1 -l-0.3(a- 1) « 
1-I-0.6A (jVladimirova et aril200fit) . Here pa and pb are 
the unburned and burned densities, respectively. Note 
that when a = 1, we have S = Sq and 1 = Iq. An em- 
pirical calibration was made based upon flames propa- 
gated with the multistage flame (described below), find- 
ing a « 1-I-1.5A. We attribute the difference to the en- 
ergy release occurring in the second stage, slightly af- 
ter the ADR flame front. The front thickness is set 
to be several grid points (spaced by Ax) per interface 
I = bAx, usually with b = 4. Further details of the imple- 
m entation of this flame-cap turing method may be found 
in lVladimirova et alJ l|2006j) . 

The above description of the ADR flame capturing 
scheme describes the method of propagating a model 
flame at a specified speed. Flames passing through stellar 
material are subject to fluid instabilities, particularly the 
Rayleigh- Taylor instability, which increases mixing and 
therefore the effective flame speed. Analyzing the bal- 
ance between generation and destructio n of flame s urface 
in a steady turbulent burning regime, iKhokhlovl ^199^ 
showed that the flame surface packed in a given vol- 
ume is inversely proportional to the laminar flame speed. 
Consequently, the total burning rate, which is equal to 
the product of the laminar flame speed and the flame 
surface area, is independent of the laminar flame speed 
in this regime. Numerical calculations of flames propa- 
gating in a vertical column demonstrate that the effec- 
tive flame speed is determined by the la rge-scale niotion 
rather than by the local flame properties l)Khokhlovll99!Tl: 
iZhang et al.ll2006|l . To account for flame acceleration on 
scales unresolved by the simulation, we use a turbulent 
subgrid flame speed l|Khokhlovlll99fit) . 

5sub = 0.5v/A^Ax, (4) 
where g is the local gravitational acceleration. As an 
input to ADR model we take 5 = max(5'iani,5'sub), where 



the laminar flame speed, 
5,a™ = 9.2 X 10^ 



0.803 

^ cms-\ (5) 



2 X IC g cm 



is from iTimmes fc Wooslevl(p9l . 

The subgrid flame speed (eq. 0]) depends on the At- 
wood number A. To compute A exactly, one must know 
the densities ahead of and behind the flame interface. 
The straightforward approach of "taking a probe" on 
both sides of the flame works well with models in which 
an infinitely-thin interface is reconstructed inside the 
computational cell and information about burned and 
unburned fluid is available in every partially burned cell. 
In the ADR flame model, to compute A in the partially 
burned cell, one must "take a probe" at the location sev- 
eral computational cells away. It is not easy to develop 
a robust algorithm for locating pure fuel and pure ash 
in the vicinity of a particular partially burned cell. Ef- 
ficiently implementing such an algorithm in a parallel, 
domain-decomposed simulation code such as flash is 
even more complicated because the information needed 
to perform the calculation for a given computational zone 
may not be directly available. 

An alternate approach is to use local information to 
estimate the burned and unburned s t ates from the par- 
tially burned state (jKhokhlovl 119951 iVladimirova et all 
l2006j) . The approach is based on the proportionality be- 
tween the mass fraction of product (f) and the fraction of 
energy released, so that the thermodynamical state can 
be parameterized by (j). In particular, the density of the 
partially burned fluid can be obtained from the unburned 
state and d>. 



P = 



Pu 



l-l-(a-l)0' 



1-1- 



7-1 gPu 
7 Pu ' 



(6) 



Here, as above, a = pu/ Ph is the density ratio across the 
interface and q is the mass-specific energy released in 
the fiame. We approximate the EOS with 7, defined 
by f = (7— \ )p£ , where £ is the specific internal energy. 
The above relations assume a very subsonic flame, but 
A is not required to be small. Under these conditions, 
the pressure variation across the flame is insignificant, 
p K Pu, and equations © can be solved for a and pu, 



Pu- 



p 



\-eo 



a=l + -^, where e=2_J.^^. (7) 

l-e<p IP 



Then the Atwood number A = (a - l)/{a + 1) and the un- 
burned density pu needed to evaluate the flame speed may 
be computed using eq. ((TJ. For small Atwood numbers. 
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Fig. 6. — Diagram of three-stage model flame front structure for a 
flame propagating to the right. The upper panel shows the progress 
variables used to track the flame activity: ipl models the '^C flame 
front, (/)2 models progress to NSQE, and <j)3 models relaxation to 
NSE. The lower panel demonstrates the elements used to represent 
the stages and their relation to the progress variables. 



we note A « e/2, which serves as a phi- independent ap- 
proximation. Similarly, 5iam is calculated using the local 
density rather than p^. This implementation has been 
tested for a variety of conditions and was found to per- 
form quite well in situations without extreme shear. 

3.2. The Stages of Nuclear Burning and Energy Release 
in the Model Flame 

As mentioned in section [21 the structure of the nuclear 
flame can be modeled as a three-stage process in which 
each consecutive stage occurs on a longer timescale. In 
our model, the progress of each of these stages is tracked 
with progress variables 4>i, 4>2, and 4>3j which evolve from 
to 1. This is shown schematically in the upper panel 
of Fig. El where the flame is propagating to the right. 
Material begins as a mixture of '^C and '^O. In stage 1, 
the '^C is consumed, in stage 2, the resulting ash burns 
to NSQE; and finally in stage 3 material relaxes to actual 
NSE. The NSQE state is reached much more quickly than 
the final NSE state, and most (> 90%) of the energy re- 
leasc in burning to NSE is released upon reaching NSQE 
llKhokhlovll 1 995j) . Having separate progress variables for 
NSQE and NSE allows the energy release timescale to be 
properly modeled while still tracking when material will 
reach actual NSE abundances. The resulting NSE state 
is not static, but can change as the material evolves hy- 
drodynamically (e.g. decompression in a rising bubble), 
and due to weak interactions. The treatment of this post- 
fiame NSE evolution is described in section 

Because the burning of '~C occurs on much less than 
the grid scale, and the physical scale of this flame front is 
an important factor in setting its propagation speed, the 



ADR flame scheme is used to evolve 01 with the param- 
eterized flame speed discussed above. This reflects the 
propagation of the actual flame front through the star. 
Although the consecutive stages separate cleanly on a 
microscopic level, as indicated in figures 1 2141 in order 
to allow for proper advection of the unresolved transi- 
tion fronts in the model, all three stages are allowed to 
be partially progressed in the same cell. However, later 
flame stages only progress when previous flame stages 
have completed. When < 1, (/)2 = 0, but behind the 
propagating '~C flame, when (j>i = l, 

• (8) 



tnsqe 



Here tnsqe is the timescale for reaching the NSQE state 
evaluated in nuclear network studies discussed in detail 
in the following section. The speed of this evolution and 
that for 03 below are limited such that 02 will change by 
at most 2% in a single explicit hydrodynamic timestep, 
in order to reduce noise. 

As noted earlier, the bulk of the energy released in the 
burn to NSE is released upon reaching NSQE. Accord- 
ingly, and for simplicity, all the energy available in the 
initial burn to NSE is released in the first two stages, so 
that the nuclear energy release in the flame is given by 

( (i>iAQc + MAQNSE-AQc) when AQnse > AQc, 

^nuc ~ \ 

[ 0iAgNSE when AQnse < Age- 

. .(.9) 

Here AQc is the energy release for burning the initial 
abundance of '^C to ^"^Mg and AQnse is the energy re- 
leased for burning to NSE. AQnse is a function of the 
local density, and is evaluated in detail in the next sec- 
tion. Even at high densities, AQnse > 3 x lO'^ergsg"' 
(see Fig. (HJ. The '^C burning releases AQc = 2.78 x 
10'''ergsg"' for Xc = 0.5, and, accordingly, the second 
form in equation (jSJ is not used in the simulations pre- 
sented in this work. 

Finally, after NSQE is reached, the final progress vari- 
able is evolved. When 02 = 1 , 

03 = i^, (10) 

tnse 

where tnse, evaluated below, is the timescale on which 
full relaxation to the NSE abundance distribution occurs. 
The evolution of the NSE state is delayed until the third 
progress variable reaches unity. Our intention is to allow 
a clear structure of the flame at low densities where the 
stages will separate significantly, and where post-flame 
NSE evolution is expected to be of relatively minor im- 
portance. At high densities the NSE timescale is short 
enough that the post-flame evolution of the NSE state 
described below in section takes over within a few grid 
cells of the ADR flame stage (0i). Note that because 
02 and 03 are not ADR variables, the flame is not sig- 
nificantly thickened by the addition of multiple burning 
stages. 

By inspecting the evolution shown by the isochoric self- 
heating calculations in section (21 we find that the abun- 
dance evolution can be approximated by selected abun- 
dant elements in the intermediate stages. The stages 
then correspond to the following transitions: stage 1, 
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NSE. This allows the abundances to be related to the 
progress variables, giving 

X,2c = (l-0i)X°, (11) 



dE _ P dp 

—J- — £nucl "I 7 -77 

dt p- dt 



(18) 



X,6o = (l-02)(1-^C), 



= (01 



(12) 
(13) 



X2,si = (02-03), (14) 
</'3 = ^^/,NSE, (15) 

where Xq is the initial abundance of '~C and the X/.nse 
consists of the elements that make up the NSE abun- 
dances. The relationship between the progress variables 
and the abundances, and their evolution, is demonstrated 
in Fig. El The distribution of material among the NSE 
elements which together make up (pj, is found from 

Xi = <i>A^{p) (16) 
and advecting the results in the usual way. Here the 
Xf^^^ip) result from the same self-heating calculation 
used to find AQnse, described in the following sec- 
tion. This also provides a reasonable crossover to the 
post-flame evolution of the NSE clement abundances de- 
scribed below. NSE material with approximately the cor- 
rect density dependent binding energy is monotonically 
produced by the flame. 

4. INPUT TO THE FLAME MODEL FROM NETWORK AND 
NSE CALCULATIONS 

The three stage model flame front (Fig. ^ requires a 
priori information characterizing the burn. In this sec- 
tion, we discuss the calibration of the energetics and 
timescales of the flame model. Using sclf-heating cal- 
culations allows us to make progress by removing the 
constraint on the size of time steps that occurs with 
an explicit hydrodynamics method (described above). 
Our principal approach for determining energetics and 
timescales is to perform self-heating simulations with the 
large nuclear reaction network, and this method works 
for most of the densities of interest. These self- heating 
calculations must still resolve the temporal evolution of 
the burn, which for the network calculations is accom- 
plished by employing a variable time step ordinary dif- 
ferential equation integrator. Although requiring small 
timesteps during some of the evolution, the calculations 
are tractable. At the lowest densities, however, the rel- 
atively long burning timescale forces us to obtain the 
abundances and binding energy of the NSE state by 
applying direct NSE calculations to the results of the 
self- heating network calculations. These calculations are 
presented in this section, and the direct NSE method 
is described in Appendix As mentioned above, the 
calibration presented here applies to the case of a 1:1 
mass ratio of C and O. The model flame may be applied 
to any constant mass ratio mixture, but would require 
the calculations described in this section to be repeated 
for the particular case (and use of the appropriate input 
flame speed). The effect would be a different calibration, 
with longer or shorter timescales and a different energy 
release. 

The network we use contains 200 nuclides up to Kr, 
shown in Fig. [3 We solve the evolution of the species 
and heat equation in time, 

— -— f-/'' -Fr® ) (17) 



where Yj is the abundance of species / in mol g~' , and E 
is the internal energy per gram. The term '"dest('"prod) 
the destruction (production) rate per volume of species i 
owing to thermonuclear reactions, which is in fact a sum 
over all the reaction paths in which the nuclide partic- 
ipates. Cnuci is the energy generation rate due to ther- 
monuclear burning, a sum of the reaction Q-values for 
each pathway times their rates per gram of material. We 
solve for the evolution of the nuclide abundances and 
heat equation, cq. H17|) and cq. H18|) . using an operator- 
split formalism. For constant pressure, eq. H18|) can be 
written as 



dT 
dt 



■ ^nucl ' 



(P_dE\dP/dT 
"'^^~~^p)WJd~p 



(19) 



where all partial derivatives are either at constant T or 
p, and Cv = If is specific heat in constant density. In 
writing cq. (|19|l . we neglect the terms of change in inter- 
nal energy due to the change in the species. These terms 
are quite small compared with enuci except right at NSE. 
Eq. H17() is stepped forward in time at a constant T and 
p using a semi-i mplicit scheme combined with a sparse 
matrix solver fsee lTimmeslll999l and references therein). 
T is updated by finite differencing H19|) and solving cx- 
phcitly. When dp/dt = (isochoric), eq. lfTB|) is simplified 
to be dE/dt = Enuci, i-e., all the energy input from nuclear 
burning is used to increase the internal energy and hence 
the temperature. In this case, we update T from the 
updated E via nuclear burning by calling the equation 
of state, which takes into account the contribution due 
to the change in the species. For constant P (isobaric), 
some part of the burning energy is used to expand the 
material. Using the constant P and updated T f eg. 119(1 . 
p can be updated by calling the equation of state. 

The thermonuclear reaction rates are taken from an ex- 
panded version (Schatz 2005, pr ivate communication) of 
the rate compilation REACLIB l|Thielemann et al.lll98(]l: 
iRauscher fc Thielemann! |2000|) . We have also included 
the temper ature-dependent nuclear partiti on functions 
provided bv IRauscher fc Thielemannl l|2fl0fD . both in the 
determination of the rates of inverse reactions and in 
our determination of NSE abundance patterns. To en- 
sure that the rates obey detailed balance in NSE, the 
inverse reaction rates are derived directly from the for- 
ward rates based upon equilibrium equations that ac- 
count for Coulomb interactions on the chemical poten- 
tial (see Appendix El-® Our calculations demonstrate 
that the average binding energy and the abundances 
of the final NSE state from the self-hcating calculation 
agrees with the results of direct NSE calculations to 
within 1%. We incorporate the effects of electron screen- 
ing of thermonuclear reaction rates, adopting the rela- 
tio ns for weak screening and strong screening provided 
by iWallace et al. (1982.) . Finally, contributions from 
weak reactions are included u s ing the rate s provided by 
iLanganke fc Martmez-Pinedol p000ll200l|) . 

^ Common practice is to use the reverse rates included in the 
REACLIB tabic, rather than calculating them explicitly from the 
forward rates. Note, however, that inverse reactions for '^C and 
'^0 burning are currently omitted from REACLIB. 
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Neutron number 



Fig. 7. — Set of nuclides included in the 200 nuclide network 
(boxes) and additional neutron rich nuclides added for the NSE 
calculations (crosses). Shaded boxes indicate stable nuclides. 



We apply this network in a collection of representative 
self-heating calculations, which are intended to inform 
our flame model both of the timcscales required to reach 
NSQE or NSE and of the accompanying nuclear energy 
release. As discussed in §12 an isochoric calculation pro- 
vides a good representation of the burning in an actual 
flame. We also use isochoric calculations below to discuss 
the importance of screening. The initial temperature is 
set as Tg = 1.7, which is sufficient for '^C and '^O ignition 
but lower than the final temperatures obtained at NSE 
for the densities of interest. As with all calculations fea- 
tured in this work, the initial composition was taken to 
be half '^C and half '^O by mass. The calculations were 
continued until NSE was achieved, except at lower den- 
sities (p < 10^ g cm"^), for which the timescales to reach 
NSE are extremely long. Although the self-heating cal- 
culations at very low densities do not yield the correct 
NSE composition and released energy, they do provide 
temperatures close to those of the NSE states. Hence, 
we apply the temperature achieved from the self-heating 
calculation to the direct NSE calculation to get the cor- 
rect released energy and composition. 

The results of these systematic studies are presented 
in Table 01 and 01 where, for the critical range of densi- 
ties, we tabulate the temperature achieved, the energy 
release in erg/g, and the mass fractions of "^He and pro- 
tons. Th e results from t wo other networks, aproxlQ and 
torch47 l|Timmesl 11991 . are listed as a comparison for 
the isochoric case. As a reference point, the conver- 
sion of an initial composition of half '^C and half '^O by 
mass (with a mean binding energy per nucleon of 7.829 
MeV/nucleon) to pure ^*'Ni (8.643 MeV/nucleon) yields 
0.814 MeV/nucleon or equivalently 7.85 x lO'^ crg/g. 
The lower values of the energy release contained in Ta- 
ble reflect the presence of significant concentrations of 
^He (7.074 MeV/nucleon) and protons; these are charac- 
teristic of NSE distributions at higher temperatures, and 
can be seen to increase with increasing density and tem- 
perature. Reliable determinations of the compositions 
under these conditions are clearly crucial to our obtain- 
ing a proper measure of the energy release, as described 
in our discussion of the implementation of our flame en- 




2 - 

10^ ID'* 10" 10'° 

Density [g cm"'] 

Fig. 8. — Energy released and NSE temperature for burning 
C/O as functions of density of the fuel for isochoric self- heating 
calculations with (solid lines) and without (dashed lines) electron 
screening and Coulomb corrections. For comparison, the contour 
corresponding to the condition (for Ye = 0.5) that the nuclear sta- 
tistical equilibrium consists of equal parts by mass of '^He and '*Ni 
is shown (dotted line). 




Density [g cm"'] 

Fig. 9. — Final (NSE) abundances of select species after C/O 
burning for the cases shown in Fig. |3 Dashed and solid lines have 
the same meaning as there. 



ergetics. 

The "ignition line", defined as the temperature-density 
relation of the final NSE state from the self-heating cal- 
culations, is presented in Fig. [HI This is shown for 
the isochoric case in order to facilitate the discussion 
of Coulomb screening corrections below. At all of the 
densities shown the electrons arc fully degenerate. For 
comparison, we sketch the contour in the logp— plane 
corresponding to the condition (for Ye = 0.5) that the 
nuclear statistical equilibrium consists of equal parts by 
mass of "^He and ^^Ni. The greatest proximity of our ig- 
nition line to the equal mass "^He/^^Ni contour is seen to 
occur at the highest densities (e.g. 10^ gem"-') and corre- 
spondingly high temperatures. The free particle densities 
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realized for the case p = 10^ are considerably lower, the 
matter is dominated by iron-peak nuclei, and the corre- 
sponding mean binding energy per nucleon most closely 
approaches that of pure ^^Ni. Note that our ignition 
line can approach - but never overlay - the equal mass 
"^He/^^Ni contour for a very straightforward reason. Re- 
call that our initial condition consisted of equal masses 
of '~C and ^^O; the mean binding energy per nucleon 
for this mixture is ^ 7.8 McV/nucleon. However, the 
mean B/A for a mixture of equal parts by mass ^He and 
^^Ni is also B/A~7.8 MeV/nucleon. It follows that the 
*He-^^Ni contour is physically not achievable for this fuel 
composition. 

4.1. Coulomb Correction and Energetics 

The impact of our inclusion of Coulomb corrections 
and screening effects both on the NSE temperatures 
achieved and on the energetics are illustrated in Fig. |H1 
A detailed discussion of the manner in which Coulomb 
and screening effects are utilized in the determination 
of screened reaction rates and nuclear statistical equilib- 
rium configurations is presented in Appendix El In or- 
der to facilitate comparison, the case without screening 
neglects both reaction rate screening and Coulomb cor- 
rections to the equation of state (EOS). This necessitates 
an isochoric comparison due to the lower pressure with 
Coulomb corrections in the EOS. We note from the figure 
that even at the highest densities, where screening effects 
are most pronounced (rse^i =Z^l^e^{'^'Kne/3)'^l^ /kT = 10), 
the differences in energy release is negligible and the dif- 
ferences in temperatures achieved are quite small. For 
the limiting case oi p = lO'" g cm~^, the temperature 
achieved without screening and Coulomb correction is 
T) = 10 while with screening it is = 10.5, increasing by 
5%. For the same conditions, the difference in energy 
release is less than 1%. 

In order to understand how plasma Coulomb correc- 
tions affect the final temperature and energy release, it 
is useful to look at how the final (NSE) composition 
changes in detail. Fig. IHl demonstrates the difference of 
the NSE compositions with screened {solid lines) and 
unscreened (dashed lines) network calculations. The dif- 
ference is again quite modest. Generally the inclusion of 
plasma Coulomb corrections energetically favors assem- 
bling fewer, higher charge ions from lower charge ones. 
Consistent with this, we find that the abundances of 
the high-charge species increase at the expense of low- 
charge ones. In light of the minor change in abundances, 
the ideal-gas contribution to the heat capacity (which 
is related to the mean molecular weight of the ion gas) 
changes very little. Thus the the difference in temper- 
ature at the high density between the screened and un- 
screened cases is likely due to the difference in the EOS 
being used. The Coulomb corrections decrease the heat 
capacity of the ion gas, allowing higher temperatures to 
be reached for the same energy release. 

The degree of agreement between the energy release 
with and without screening terms is remarkable. A clue 
to the origin of this match is that there is one exception to 
the general rule that Coulomb enhances heavy element 
abundances: protons. Although the ^He, '^C, and '^O 
abundances are decreased in the screened case, protons 
actually increase. This is because we are working with 
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Fig. 10. — Nuclear burning timcscales versus ash temperature 
of C/O fuel calculated by 200 nuclide self- heating network. Shown 
are curves for the timescales to NSQE and NSE for isochoric {solid 
lines) and for isobaric (dashed lines) one-zone calculations. 



the constraint that = 0.5, while Ni is by far not the 
most abundant Fc-group nucleus, with more neutron-rich 
nuclides like ^^Co and ^^Ni being favored. Thus it ap- 
pears that although high-charge nuclei are favored with 
screening included, this is offset in energy release by the 
concomitant increase in proton abundance. 

4.2. NSQE and NSE timescale calculations 

As noted above, the three-stage flame model requires 
timescales for the second and third stages of burning. 
These timescales must be determined in advance of sim- 
ulations with the three-stage flame model, and accord- 
ingly we investigated determining these timescales with 
the self- heating networks available to us, and the effect 
of screening on these timcscales. 

The timescale to NSQE was defined as the elapsed 
time in a simulation between the point at which 90% 
of the C has burned and the maximum abundance of 
^^Si. Defining the timescale to NSE requires a little more 
care, and we investigated several approaches. In all cases, 
the timescale was measured as the elapsed time from the 
maximum abundance of ^'^Si. The results of these stud- 
ies are shown in Table |21 In Table |21 wc list the NSE 
timescales defined in several ways: the time when 90% 
of the maximum ^^Ni abundance is reached; the time 
when 90% of the maximum abundance of the most abun- 
dant Fe-peak isotope is reached; the time when 90% of 
the maximum abundance of NSE nuclei is reached. NSE 
nuclei include ''He and isotopes of Cr, Fe, Co, Ni, and 
Zn. At the lowest density end, Ni dominates the NSE 
compositions, so all the approaches yield the same NSE 
timescales. At the highest density end, the dominant 
NSE composition is a mixture of He and Fe-peak iso- 
topes. By including Fe-peak nuclides other than Ni, we 
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TABLE 2 

NSQE AND NSE TIMESCALES FROM SELF-HEATING NETWORK CALCULATIONS WITH THE 200 NUCLIDE NETWORK. 



max. 2^Si abund. 90% of ^^Ni abund. 90% of Fe-peak nucleus 90% of NSE nuclei 
log(p) '^NSE "HMSQE TNSE TNSE "T^SE 

(gcm-3) (10" K) (s) (s) (s) (s) 



7.0 


4.56 


2.21E-03 


2.90E+00 


2.90E+00 


2.86E+00 


7.1 


4.75 


5.46E-04 


6.62E-01 


6.62E-01 


6.65E-01 


7.2 


4.95 


1.39E-04 


1.70E-01 


1.70E-01 


1.68E-01 


7.3 


5.13 


3.75E-05 


4.94E-02 


4.94E-02 


4.56E-02 


7.4 


5.32 


1.07E-05 


1.19E-02 


1.19E-02 


9.40E-03 


7.5 


5.50 


3.28E-06 


4.33E-03 


7.16E-04 


3.02E-03 


7.6 


5.68 


1.08E-06 


1.72E-03 


2.26E-04 


1.02E-03 


7.7 


5.85 


3.81E-07 


7.80E-04 


8.08E-05 


3.85E-04 


7.8 


6.01 


1.45E-07 


3.61E-04 


3.22E-05 


1.46E-04 


7.9 


6.17 


5.85E-08 


1.68E-04 


1.45E-05 


5.64E-05 


8.0 


6.33 


2.50E-08 


8.01E-05 


7.08E-06 


2.27E-05 


8.1 


6.50 


1.12E-08 


3.87E-05 


3.66E-06 


9.56E-06 


8.2 


6.66 


5.19E-09 


1.89E-05 


1.97E-06 


4.22E-06 


8.3 


6.82 


2.50E-09 


9.41E-06 


l.lOE-06 


1.94E-06 


8.4 


6.99 


1.26E-09 


4.77E-06 


6.01E-07 


9.35E-07 


8.5 


7.15 


6.58E-10 


2.48E-06 


3.45E-07 


4.67E-07 


8.6 


7.33 


3.58E-10 


1.30E-06 


2.03E-07 


2.41E-07 


8.7 


7.50 


2.02E-10 


7.00E-07 


1.20E-07 


1.28E-07 


8.8 


7.68 


1.19E-10 


3.82E-07 


7.17E-08 


7.01E-08 


8.9 


7.87 


7.13E-11 


2.12E-07 


4.32E-08 


3.91E-08 


9.0 


8.06 


4.42E-11 


1.21E-07 


2.52E-08 


2.24E-08 
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Note. — The NSQE timcscalc is defined as tlie elapsed time in a simulation between the point at which 90% of the 
C has burned to the maximum of "^Si abundance. The time scales to NSE were determined as the elapsed time from 
the maximum abundance of ^^Si, to the 90% of the maximum ^^Ni abundance, to the 90% of the maximum abundance 
of the most abundant Fe-pealc nucleus, and to 90% of the maximum number of NSE nuclei, respectively. NSE nuclei 
were '^He and isotopes of Cr, Fe, Co, Ni, and Zn. 



can see that the timescales differ at most within a factor 
of 3 at the high density end. Inchiding He decreases the 
timcscalc further by a factor of 1.5 at the high density 
end. In this paper, we choose the last definition as NSE 
timcscalc in implementing the flame model. The NSQE 
and NSE timescales from the self-hcating calculations arc 
fitted by 

7-NSQE = exp(182.06/r9-46.054) (20) 

and 

TNSE = exp(196.02/r9-41.645) . (21) 

We express the timescales as functions of final temper- 
ature, which is determined by the final density and en- 
ergy release in self-heating calculations. This is a quite 
robust parameterization as seen from the comparison of 
timescales determined from isochoric and isobaric cal- 
culations in Fig. Uni We note that our N SQE/NSE 
timesca les are consistently lower than those of lKhokhlovl 
l)1991a() . but we are unable to determine the reason for 
this difference. 

The influence of screening on the burning timescales, 
as inferred from our self-hcating calculations, is revealed 
in Fig. ^2 Once again we use an isochoric case to isolate 
the effects of rate screening. In contrast to our find- 
ings for temperature and energy, the impact of the in- 



clusion of Coulomb and screening effects is seen here 
to decrease the timescales to achieve NSQE and NSE 
by factors of 2-3. This reflects the level of individual 
screening enhancements for the critical reactions that 
govern the flows into the iron-peak region. An isobaric 
burn, representative of what will occur in the star, is 
also shown, and has much longer timescales at lower 
densities due to the lower temperature on expansion. 
For comparison, the dynamical timescale for homolo- 
gous collapse, Tdyn = (247rGpr'/2 = 0.14(10^ g cmr"' / pfl^ s, 
is superposed [dotted line) on the figure. Note that 
tnse > Tdyn for /9 < 4 X 10^ g cm"^ and tnsqe > Tdyn for 
p < 1.5 X 10^ g cm"\ 

5. POST-FLAME ENERGY RELEASE AND 
NEUTRONIZATION 

Once the fiame has passed through a piece of stellar 
material, the remaining ashes have burned to NSE, which 
is predominantly comprised of iron-peak nuclei, "'He, and 
free protons. This NSE composition, however, should not 
be considered static. It changes as the density decreases, 
either in a rising bubble or due to the expansion of the 
star, resulting in a varying ionic heat capacity and more 
importantly causing the difference in nuclear binding en- 
ergy between the initial '^C/'*0 fuel and the NSE state 
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Fig. 11. — Nuclear burning timcscalcs versus density of C/O fuel 
calculated by 200 nuclide self-heating network. Shown arc curves 
for the timcscalcs to reach NSQE and NSE with {solid lines) or 
without {dashed lines) screening and Coulomb corrections for an 
isochoric burn and with corrections for an isobaric burn {dash-dot 
lines). Superposed is the dynamical timescale {dotted line), r^y^. 



to vary during the course of the hydrodynamic evohition 
of the supernova. During this extended evolution, there 
is time for weak interactions to take place in the ash 
material. Dominated by electron captures, these pro- 
cesses lower F(, and emit neutrinos. Below we discuss 
these two physically distinct effects which bring about 
changes in the NSE abundances. The method used to 
calculate (Coulomb-corrected) NSE is described in Ap- 
pendix ^ with the set of nuclides used in our network 
calculation supplemented with neutron-rich Fe peak ele- 
ments as shown in Fig. 

5.1. NSE adjustment to the evolving state variables p 

and T 

The hydrodynamic evolution of the ashes, i.e. the ex- 
pansion following the burn, causes a decrease in p and 
r on a hydrodynamic timescale, which in turn causes a 
change in the NSE mass fractions as the nucleons reor- 
ganize to maximize entropy for the new thermodynamic 
state. As long as the hydrodynamic timescale substan- 
tially exceeds the charged-particle (strong) nuclear reac- 
tion timescale, the approximation of instantaneous read- 
justment of the composition of the ashes to the NSE ap- 
propriate for the local thermodynamic state is a good 
one. 

The electron fraction, defined as Y„ = V Z, — / V A, — , 

where m,- is the nuclear mass of species of the '^C/'^O 
mixture is 0.5. Weak interactions are assumed to be 
negligible during the burning phase of the flame since 
the weak interaction timescale far exceeds the burning 
timescale. At first, the composition of the hot, unex- 
panded ashes has a significant fraction of the mass in 
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Fig. 12. — Contours of binding energy per nucleon, in MeV 
per nucleon, for the NSE state with Y/, — 0.5. Also shown is as a 
function of density the final temperature reached (dot-dashed line) 
by the 200 nuclide self-heating network calculation starting from 
an equal mixture of C and O by mass and burning to NSE. The 
line of 50% '*Ni mass fraction (dotted line) in a two-component 
NSE C^He and ^^Ni) indicates the temperature above which light 
particles dominate. 



free protons (0 MeV) and "^He (7.074 MeV/nucleon). 
For example at = l,Xp = 0.024, Xhb = 0.095 and at 
p% = 10, Xp = 0.027, Xhc = 0.135. The remainder predom- 
inantly consists of a mix of strongly bound Fe-peak nu- 
clei, dominated by ^^Fe. The nuclear binding energies 
reached in the hot ashes immediately following the burn 
are therefore significantly smaller than the binding ener- 
gies obtained by a cold NSE dominated by ^^Ni and lie 
along the dot-dashed line in Fig. [T^ 

The binding energy gain per nucleon of such a mixture 
compared to the initial '-C/'^O fuel, which has a binding 
energy per nucleon of 7.828 MeV, varies with density 
but remains for > 1 in the vicinity of 0.4 MeV. As 
the material expands and cools, NSE progressively favors 
heavier nuclei. The reason for the shift from an alpha rich 
NSE to an alpha poor NSE during an expansion is that 
lower entropies allow for a reduction of the total number 
of particles, since the number of accessible microstates 
of the ensemble has decreased. The light particles are 
reassembled into fewer but more tightly bound Fe-peak 
nuclei. At constant Ye = 0.5, the composition would reach 
virtually 100% ^^Ni at freezeout temperatures. In other 
words, at constant 7^. = 0.5, the binding energy would shift 
from ~ 8.2 to 8.643 MeV per nucleon and an additional ~ 
0.4 MeV of binding energy would be released per nucleon 
during the expansion phase. This is a significant fraction 
of the total energy released and must be accounted for. 

Furthermore, because the average nucleon number per 
nucleus A = y'. A/ — / V,. — grows as the NSE distribution 
is further shifted with decreasing temperature toward the 
iron-peak nuclei, the ionic heat capacity decreases. Be- 
cause we have conditions that are at the borderline of 
electron degeneracy, we estimate this to have a small ef- 
fect. 

5.2. Neutronization 

Virtually all extant one-dimensional calculations of 
SNe la explosions find that the innermost regions (the in- 
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Fig. 13. — Contours of time derivative of electron fraction, Ye, 
due to electron and positron cap ture and emission, evaluated at 
Ye = 0.5. Other lines are as in Fig. 1121 
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Fig. 14. — Contours of time derivative of electron fraction. Ye, 
due to electron and positron capt ure and emission, evaluated at 
Ye = 0.48. Other lines are as in Fig. ^ 



ner « 0.2 Mq) suffer sufficient changes in Ye, due to elec- 
tron captures, to sig n ificantly change t h e compositions 
llNomoto et all 119841: iKhokhlovl llQQlbl: llwamoto etHl 
Il999j) . The use of recent electron capture rates, based 
on shell model diagonalization, reduces the neutron ex- 
cess in the central regions somewhat ()Brachwitz et al.l 
l200nj) . the gross features however remain. Neutroniza- 
tion is a consequence of the high Fermi energies and the 
correspondingly high rates of electron capture reactions 
under these conditions. The time rate of change of the 
electron fraction Y^, also referred to as the neutronization 
rate, and the amount of energy carried away by neutrinos 
per gram of stellar material per second E'^ are calculated 
as NSE av erages by convolving tabulated wea k interac- 
tion rates l)Langanke & Martmez-Pineddl200ll) with the 
NSE abundances. We consider electron decay, positron 
captures, electron captures and positron decay of pf-Shell 
nuclei as well as free nucleons, so that 

Ye=y^(xf+xr-xr->^f) (22) 

^^-^ mi \ / 

E-=y^(xf+xr) (23) 

^-^ nv, \ J 

where we are following t he notation from 
iLanganke fc Martmez-Pinedcl l)200lD in which 
Xed ^XP' ,\'' ,\P'' are given in units of s"' and ,Y 
are given in units of MeV s"' . Neutrinos emitted are 
taken to be free-streaming. 

At high densities (p9 ^ 1) the NSE averaged neutron- 
ization rate for self conjugate matter is rather fast with 
-0.05 s-i > Ye>-0.5 s-i (see Fig. US)). The rate of neu- 
tronization is strongly dependent on the electron fraction 
(compare Figs. 1131 and lT^ . since electron capture rates 
are overall decreasing with increasing neutron number for 
a given element. Due to the strong decrease of the rate 
of neutronization with decreasing electron fraction, one 
typically finds the composition of these inner regions to 
consist of moderately neutron-rich, strongly bound iron- 
peak nuclei, viz; ^^Fe, ^"^Fe, and ^^Ni. 

If neutronization occurs rapidly enough to lower Y^ to 
~ 0.48, then the binding energy per nucleon of the neu- 
tron rich NSE state compared to the NSE state with 
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Fig. 15. — Contours of binding energy per nucleon, q, in MeV 
per nucl eon, for the NSE state with Ye = 0.48. Other lines are as in 

Fig. im 



Yc = 0.5 is increased by ~ (0.2 - 0.3) MeV for a range 
of given temperatures and densities (compare Figs. 1121 
and US)). In the cold NSE limit, where ^^Pe (8.736 
MeV/nucleon) and ^'^Ni (8.732 MeV/nucleon) are the 
dominant nuclei of the Y^ = 0.48 NSE abundance distribu- 
tion, the final binding energy per nucleon increase com- 
pared to NSE of self conjugate matter, which is entirely 
dominated by ^^Ni, is on the order of 0.14 MeV. Not all 
of this energy, however, remains available to expand the 
star. Neutrino losses accompanying neutronization im- 
pact the energetics: some fraction of the energy release 
(being somewhat dependent upon the density and corre- 
sponding electron Fermi energy) will be lost in neutrinos 
(see Figs. [Tlfc[T7|). 

Another important effect of neutronization lies in its ef- 
fect on the Fermi level of the electrons. As electrons get 

captured, Yg decreases and since relativistic electron de- 

4 

generacy pressure scales as (pFc)^ this results in a change 
in the ratio of electron degeneracy pressure to ionic pres- 
sure. Because we are close to the borderline of degener- 
acy, it is difficult to estimate the impact of this effect on 
the thermal histories of tracer particles being followed in 
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Fig. 16. — Contours of neutrino energy loss rates, e^, during 
electron and positron cap ture s and emissions, evaluated at Yg =0.5. 
Other lines are as in Fig. 1121 
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Fig. 17. — Contours of neutrino energy loss rates, e^, during 
electron and positron capt ure s and emissions, evaluated at Ye = 0.48. 
Other lines are as in Fig. 1121 



the simulation. 

The impact ol neutronization on the rise times ol hot 
bubbles and the subsequent growth of the Rayleigh Tay- 
lor instabilities could help to distinguish between com- 
peting ignition scenarios. In addition, it affects the re- 
lation between the composition of the progenitor white 
dwarf and the mass of ^''Ni synthesized in the explosion 
<|Timmes et al.ll2003|) . 

5.3. Reduced Set of Nuclei for Evolution in NSE 

While the NSE state is calculated using many nu- 
clear species, it is not feasible to store the abundance 
of each of these at every grid point. In reality only a 
small number of quantities which are derivable from the 
set of abundances. {X/}, arc essential to accurate hydro- 
dynamic and energetic simulation of the NSE material. 
These are the electron fraction, Yg, the average number 
of nucleons per Jon, A, the average binding energy per 
nucleon, q = Q/A, and the coulomb coupling parameter 

T = 7Fl^e^{ATm,/3fl^/kT, where rig is the electron num- 
ber density. The quantities Y^, A, and T arc required to 



determine the local pressure from density and internal 
energy content, while q provides a measure of the energy 
stored as rest mass energy of (comparatively less bound) 
light nuclei; which is released as the NSE state evolves. 

Working in the approximation that m,- = A,- a.m.u. re- 
duces several quantities to simple sums over the abun- 
dance set, Ye = IM = X^/^'M'' ^^^"^ ^ - 

'Y^-XjQi/Ai. This approximation is accurate at the level 
of 10"^ for these averages, but is unsatisfactory for the 
actual calculation of NSE due to its dependence on small 
relative differences between species. With this simple 
linear form, we see that it is possible to select a set of 
abundances {^/} which give the same Y^, A and q, and 
give the correct mixing evolution of these quantities (es- 
sential for an Eulerian code), but which consists of only a 
select few elements. This representative set can then be 
used in the simulation without losing information about 
the energetic evolution of the material, giving appropri- 
ate Lagrangian temperature histories which can be used 
for post-processing. 

Choosing this reduced set is much like finding basis 
vectors which span the 3D space {Ye,\/A,q), except that 
is constrained to be 1, so that the set must both 
span and bound the available space. We restrict ourselves 
to elements which exist in nature and thus have known 
{Z/A,\/A,Q/A). The goal of this method is to utilize 
existing hydrodynamic tools which are known to treat 
the Xi properly. An alternative is for the basic quantities 
{Ye,l/A,q) to be treat ed directly by th e hydrodynamics 
method, as is done by lKhokhlovl IpOSj) . 

The quantities Y^, A, and q are three of the f our quanti- 
ties tr acked by the nuclear kinetic equations of lKhokhlovl 
ljl995j) (and related references). The fourth is a variable 
to track the relaxation from NSQE to NSE, which they 
approximate to have no effect on the energy release, and 
which corresponds to our third flame stage. .Khokhlovl 
lj2000(l makes the additional approximation that A (their 
Fi) is constant in the ashes. 

We do not treat F in any special way. It is expected, 
however, that because the representative nuclei do not 
have very different Z from the real ones, the error is not 
too bad. NSE is generally a mixture of ^He and many 
species near ^^Ni, making T fairly well matched if all the 
other quantities are matched. Thus the benefit of devel- 
oping a 4D matching method that includes Z^/^ does not 
appear necessary for the la supernova problem, where F 
is moderate (~ 0.5), and declining, during the evolution. 

The simple dependence on a F which contains Z^/-' is also 
not a general feature of Coulomb corrections, but only of 
a particular approximation used for EOS corrections. 

The set of elements we have chosen, which includes 
several necessary for evolution of the flame, consists of 
neutrons, 'H, ^He, ^^Ne, ^^Mg, ^''Ni, ^^Ni, ^''Fe, ^^Ni, 
^"^Ni, and ^*Kr. (The full set in the simulation also in- 
cludes '^C, '^O, and ^^Si, used in the burning stages.) 
These are shown in Fig. ^] along with the full set used 
in the nuclear network and NSE calculation. The set of 
neutron-rich heavy elements is used to follow in detail the 
energy release as the material neutronizes. At most four 
species are used to construct a given {Ye,l/A,q), where 
the species are chosen by comparing these coordinates 
to the face planes of the four-sided polyhedra formed by 
four-element sets, starting from (^He, ^"^Mg, ^^Ni, ^^Ni) 
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Fig. 18. — Binding energy per nucleon QjA vs. Z/A (respectively 
ci and Ff, of a pure sample) for nuclei used in simulations. Shown are 
the 200 nuclei in the dynamic network calculation (+), and those 
that were added to extend the Ye range of the NSE calculation 
(x). Dashed lines extend toward neutrons on the left and 'H on 
the right; ^He is also off the plot. Boxes indicate those included 
in the representative set. The set of heavy nuclei is chosen to 
allow spanning of the entire possible area of q and Ye (outlined) 
and also allows accurate representation of A, precise matching of 
A is achieved in most cases in all but the lowest Ye portion of the 
solid-outlined area. 



and walking to more neutron rich elements as necessary. 
If the desired combination lies outside of polyhedra avail- 
able with these elements, A is not matched and the three 
species which form the triangle enclosing the and q of 
the state are used. The set above allows us to match Ye 
and q exactly to that of the NSE state for Ye > 0.42. A 
is also matched exactly in most cases, except for Ye near 
0.42, where it may vary from the NSE value by up to 
10%. 

5.4. Energy Evolution in NSE 

The NSE properties of the material are used to evolve 
the mass-specific thermal energy £, much like a nuclear 
network would evolve the abundances and £ between 
hydrodynamic steps. The resulting changes are either 
driven by the hydrodynamic motion, and thus occur on 
the same timescale, or are due to neutronization, which 
occurs on a slow timescale compared to the hydrodynam- 
ics. After a hydrodynamic step, the values of £„, P, T„ 
and {X/}„ are available at each grid point (« labels the 
timestep). The qn and Ye„ are calculated from {^/}„, ac- 
counting for how mixing changes these quantities. This 
state is slightly out of NSE as a consequence of the hy- 
drodynamic evolution and a time interval Af has passed 
during which some neutronization and neutrino loss has 
occurred. We calculate % and Cj, from tables using p„, 
Tn, and Ye^,,, and then set F(.,„+i =Ye.„ + YeAt. These func- 
tions are not expected to vary steeply in these variables 
compared to the amount they will change in a single hy- 
drodynamical step. We then find 7;,+i by solving 

£[Tr,+l,Ye,n+l,A(T„+i,Ye^n+l)] - q(Tn+i ,Yej,+l) = 
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Fig. 19. — Abundances from a simulation of a model flame 
passing through a 50%/50% by mass mixture of C and O at an 
initial density of 3 X lO' g cm"^. The left panel presents the abun- 
dances of selected surrogate nuclei on the entire simulation domain 
at t = 0.6s. Also shown in the left panel is the electron fraction. 
The right panel shows the detail of the thick flame with another 
set of surrogate nuclei. 

£„ - q„ + At [-YeNAinie -m„ + mp)c^ -e^], (24) 

at fixed p. Both A and q are tabulated from the NSE 
calculation using our full set of nuclear species. Finally 
this temperature is used to calculate f „+i , qn+i and A„+i , 
and the procedure described above is used to convert 
Ye.n+i; A„+i, and q„+i to {X/}„+i. 

6. ONE-DIMENSIONAL FLAME MODEL RESULTS 

In order to exercise the energy deposition and neu- 
tronization methods described, we performed a series of 
one-dimensional model flame simulations at a variety of 
densities. In each of these a flame was propagated from 
left to right through fuel consisting of equal parts by 
mass of '^C and '^O. In these simulations the input 
flame speeds were laminar flame speeds (eq. |5]); tur- 
bulent flame speeds (eq. 0) were not used. A reflecting 
boundary condition was used on the left, where the flame 
is ignited, and outflow on the right, so that fuel is free 
to move off the grid as the volume of (less dense) ash 
increases. A flame is started by setting 0i = 1 in a small 
region (with a smoothed edge) adjacent to the left re- 
flecting boundary and depositing the appropriate energy 
locally. The Figs. ll9H22l show profiles of the abundances 
for some of the surrogate nuclei and Ye at a time when 
the flame has propagated most of the way across the do- 
main. Right panels show a detail of the model flame 
region, which by definition is only several computational 
zones wide. 

Fig.s I19H21I shows the results of a model flames prop- 
agating across a lO^cm domain with fuel at densities of 
3 X 10^, 10^, and 10*^ g cm"-', respectively. The elapsed 
times of the simulations are t = 0.6, 1.2 and 6.0 s, re- 
spectively, and the simulations were performed on a uni- 
form mesh of 1280 points. The material furthest to the 
left has been burned for the longest time, and therefore 
has undergone the most neutronization. Moving to the 
right, ash is younger closer to the propagating flame. The 
surrogate nuclear set works as expected, with "^He and 
consecutively more neutron rich Fe-peak elements rep- 
resenting the NSE state. The detailed flame structure 
is also as expected, with the ''C burning being followed 
very closely (within a few zones) by conversion to ^'^Si 
representing NSQE and that on to NSE. The density 
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Fig. 20. — Abundances from a simulation of a model flame 
passing through a 50%/50% by mass mixture of C and O at an 
initial density of lO' g cm"'. The left panel presents the abundances 
of selected surrogate nuclei on the entire simulation domain at f = 
1.2s. Also shown in the left panel is the electron fraction. The 
right panel shows the detail of the thick flame with another set of 
surrogate nuclei. 
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Fig. 22. — Abundances from a simulation of a model flame 
passing through an equal mass mixture of C and O at an initial 
density of 5.0 X 10^ g cm"'. The left panel presents the abundances 
of selected surrogate nuclei on the entire simulation domain at ; = 
16.0s Also shown in the left panel is the electron fraction. The 
right panel shows the detail of the thick flame with another set of 
surrogate nuclei. 
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Fig. 21. — Abundances from a simulation of a model flame 
passing through an equal mass mixture of C and O at an initial 
density of 1.0 X 10^ g cm"'. The left panel presents the abundances 
of selected surrogate nuclei on the entire simulation domain at / = 
6.0s Also shown in the left panel is the electron fraction. The 
right panel shows the detail of the thick flame with another set of 
surrogate nuclei. 



dependence of the weak reactions is quite evident. At 
3 X 10^ g cm"-', Ye = 0.44 is reached in only 0.6 s, though 
at 10^ g cm"^, even after 1.2 s Yg has only fallen to about 
0.47. We performed an additional check on the neutron- 
ization by running a sclf-hcating calculation using the 
200 nuclide network with weak reactions. Comparing 
the degree of ncutronization, 0.5 — Fg, we found that the 
self- heating calculation at 10^ g cm"^ differed from the 
leftmost zone of a flame model simulation at initial den- 
sity of 1.2 X 10^ g cm~^ (resulting in an ash density of 
lO'' g cm-3) by 2.6% and 0.7% at 0.6 and 1.2 s respec- 
tively. Finally, the speed of the flame in the computa- 
tional cell due to the creation of the expanded ash is 
related to the laminar flame speed by v = 5iaminai;0fuei/Pash, 
and the test simulations reproduce this behavior satis- 
factorily. 

Fig. 1221 shows the result of a model flame propagating 
across a larger 8 x lO^cm domain with fuel at a density 
of 5.0 X 10^ g cm"^. As with the high density simula- 
tions, the simulation was performed on a uniform mesh of 
1280 points. The slower ncutronization time scales with 
respect to burning necessitated using the larger simula- 
tion domain to allow evolution of the flame for a longer 



time, 16.0 s. In this case, the ncutronization following 
the flame is similar to the higher density flames although 
the process of ncutronization is much slower. The de- 
tail of the flame in this case shows a broader ^^Si region, 
covering more than six computational zones, indicating 
that at this low density the flame model is beginning to 
resolve the NSE timescale. 

In addition to these results, we performed a variety 
of test simulations including model flames propagating 
away from an isolated ignition point in the middle of 
the domain. A similar calculation in which the fuel was 
moving at a constant speed across the domain gave iden- 
tical abundance profiles. The results of these tests are 
shown in Figs. |2S1 and |31 In both cases, the high de- 
gree of symmetry indicates that the direction in which 
the flame propagates from the ignition point does not 
affect the solution. The simulation presented in Fig. |^ 
tested the ability of the code to correctly propagate the 
flame in both directions as fuel moves across the grid and 
indicates a satisfactory result. 

7. CONCLUSIONS 

We have extended the ADR flame model 
ijVladimirova et all I2006|) with multiple burning stages 
and a reactive NSE post-flame state appropriate for 
simulations of C/0 deflagrations. New calculations 
were performed to calibrate flame structure, progress 
timescales, NSE properties and ncutronization. This 
calibration involved detailed calculations of several 
types, including direct flame simulation with a small 
nuclear network, one-zone isochoric self-heating nuclear 
network calculations, and direct NSE calculations. We 
have described in detail how these calculation justify 
a parameterized flame model, and what features and 
timescales must be reproduced in a large-scale simula- 
tion. By including Coulomb equation of state corrections 
to detailed balance, we established self-consistency be- 
tween screened reaction network calculations and direct 
NSE calculations. Our ncutronization is calculated with 
up-to -date weak rates |La.ngankc & Martinez-Fined^ 
12001(1 convolved with our full NSE distribution. We also 
developed a method of representing the precise binding 
energy and Y^ of the NSE state with a small set of nuclei 
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4.0 

Distance (10*' cm) 

Fig. 23. — Abundances and electron fraction from a simulation 
of a model flame passing through an equal mass mixture of C and 
O at an initial density of 1.0 X lO' g cm"^^. In this case, the "match 
head" was placed at the center of the domain so that the flame 
propagated t o th e left and right. The results may be compared to 
those of Fig. 1201 and the obvious symmetry indicates the flame is 
correctly propagating in both directions. 




6.0 

Distance (10*' cm) 



10.0 



Fig. 24. — Abundances and electron fraction from a simulation of 
a model flame passing through an equal mass mixture of C and O at 
an initial density of 1.0 X lO' g cm"^^. In this case, the "match head" 
was placed at the center of the domain so that the flame propagated 
to the left and right while the material on the grid slowl y flo wed 
to the right. The results may be compared to those of Fig. 1231 the 
principal difference being that the symmetry is centered about x = 
6.6 X lO' cm. The obvious symmetry indicates the flame is correctly 
propagating in both directions and that the slow propagation of the 
material across the simulation grid does not affect the results. 



that can then be easily treated in a general multi-fluid 
code. 

These improvements are necessary for several reasons. 
First, full-star models with embedded tracers for nucle- 
osynthetic post-processing are now feasible. Following 
the passage of the flame, the NSE state evolves as large- 
scale fluid motions change its temperature and density. 
At higher densities, electron captures make the NSE state 
more neutron-rich and lower the pressure from degen- 
erate electrons. Modeling both of these effects is criti- 
cal for the tracers to capture a high-fidelity recording of 
the temperature-density history of the burn for subse- 
quent nucleosynthesis calculations. Having a description 
of the NSE state that follows the post-explosion expan- 
sion and freeze-out is also important for evolving models 
far enough to compute lightcurves for contact with ob- 
servations. 

A second motivation for improving the flame model is 
to enable simulations of high-central density ignitions. 
Under these conditions, electron captures can play an 
important role in the energetics and yields of the ex- 
plosion, and their influence must be accurately treated 
in the model flame. Observations are increasingly find- 
ing a diversity of Type la energetics and perhaps a 
diversity of progenitors ijScannapieco fc BildstenI 120051 
^lannucci et al. 20Q^. Connecting to observations de- 
mands that simulations of the explosion also operate 
under a diverse range of conditions. With this flame 
model, we can now simulate explosions from a wider 
range of progenitors. While this paper presents results 
for a '^C/'^O mixture with a 1:1 mass ratio, the flame 
model is easily generalized to arbitrary '^C:'^0 ratios by 
varying AQc and (see eq. [H] and following discus- 
sion). The only limitation here is that Xq must be uni- 
form throughout the white dwarf. Within this constraint 
of a spatially uniform initial composition, our general- 
ized flame model is suitable for studies of thermonuclear 
deflagrations in white dwarfs under a range of ignition 
conditions. 

A study of the hydrodynamic character of precisely 
how the three-stage burning is integrated with the ADR 
flame and alternative ADR-type schemes is well under- 
way and, due to its relative complexity, will be published 
as a separate work. Here we have only presented 1- 
dimensional test calculations based on a simple integra- 
tion with the AD R flame model a s it has been presented 
in the literature ijKhokhlovll 1 995() . Future work will ad- 
dress the impact of many 3-dimensional flow character- 
istics such as flame front curvature, acoustic behavior, 
flame front stability, and the effects of flnite resolution. 
Such a hydrodynamical study is essential in order to un- 
derstand WD deflagration simulations performed with 
this flame-capturing technique, and to generally under- 
stand the behavior of flames with a reactive ash such as 
this. 

While improving the flame model is the motivation for 
this paper, our studies of explosive burning also clar- 
ify the importance of including Coulomb corrections to 
the equation of state and screening self-consistently. Al- 
though these inclusions do not signiflcantly affect the 
energy release and flnal temperature, they do decrease 
the timescales for reaching NSQE and NSE by factors of 
2-3, with the most pronounced differences at high den- 
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sities. In the physical conditions under consideration, 
these Coulomb effects are small but no n- negligible, mak- 
ing it worth evaluating them accurately. Since the effec- 
tive turbulent flame speed is independent of the burning 
rate, it is less critical to compute the local flame speed 
accurately. We caution, however, that the determination 
of the effective flame speed is based on an assumption 
of steady turbulent burning. It is therefore important to 
capture the underlying nuclear physics with fidelity. 
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APPENDIX 



PLASMA COULOMB CORRECTIONS FOR NUCLEAR STATISTICS 

At the high densities that occur near the center in the initial phases of a Type la Supernova, the screening of charged 
particle reactions by the electron background can significantly enhance the reaction rates. This is, in fact, the same limit 
in which Coulomb^ corrections to the equation of state (EOS) become important for the statistical factors which enter 
both a statistical equilibrium calculation and in determination of reverse ratios for charged particle nuclear reactions. In 
this appendix we outline our treatment of plasma Coulomb corrections which preserves consistency between screening- 
enhanced reaction network calculations used to obtain burning timescales and direct nuclear statistical equilibrium 
(NSE) calculations used to follow the later evolution of the NSE state. 

The nuclear statistical equilibrium (NSE) state is defined by /i,- = Z,/ip-l-(A, — Z,)^„ where fXi is the chemical potential 
(including the nuclear rest mass energy) of a given nuclide made up of Z,- protons and A, — Z,- neutrons, fip and /i„ 
are the chemical potentials of the free protons and neutrons. We calculate the nuclear rest mass energy as m,c^ = 
ZiiripC^ +{Ai-Zi)mnC^ -Qi, where Qi is the nuclear binding energy. To obtain the correct NSE state in a dense plasma, 
must include Coulomb corrections. The importance of these effects is measured by the Coulomb coupling parameter 
r = 7?e^ jakT where Z is the ion charge and a = (47rw/3)~'^^ is the av erage di stan c e between ions, with n rep rescntin, 
the ion number density. It has been found (|Hansen. Torrie. fc VieillefossGil977t iDewitt. Slatterv. fc Chabrienil99' 
that, to very good approximation, the thermodynamic properties of multi-component plasmas can be computed using 
the linear mixing rule, in which the correction to the extensive Helmholtz free energy is summed for the constituents 
in proportion to their numb er. This results in a correction to fij = nijC^ + ij,f + nf for which we use the fitting form 
^Chabrier fc PotekhinllT99a 



kT 



Vr,(A2+r,)-A2in( ^/^ + yi+£ 



r,-arctan a/F,- 



(Al) 



where F,- = Z; ' Fp is the ion-specific coulomb coupling parameter with Fg = e^(47r«£,/3)'''-'/A:r, A] = -0.9052, A2 = 0.6322, 
and A3 =-\/3/2-Ai/V3I. This for m is accurate for both weakly (F; < 0.1) and strongly (1 < F,- < 160) coupled liquids, 
an important feature due to the wide v arietv of Z present in an NSE calculation. We have also included temperature 
dependent nuclear partition functions ijRauscher fc Thielemannll2000t) in uf. Mor e complicated prescriptions for the 
multi-component plasma free energy have been explored ijNac^m^^fc Yudiri2005j) . but the linear mixing rule appears 
sufficient for our purposes of accomplishing correct energetics and hydrodynamics at moderate F. The abundances in 
NSE can now be found directly by using the equality of chemical potentials (detailed balance) to write Xj in terms of 
the ideal part of the free nucleons' chemical potentials, 



Xi-- 



m, 
P ' 



3/2 



exp 



z,Ai;d+/v,Xd+a-/if+z,-^c 



kT 



To solve the system, we substitute HA2|I into the constraint equations 



^X,-1 = 
V \iYe - 1 )Zi + YMi - Zi)] — = 



(A2) 

(A3) 
(A4) 



^ Note that "screening" is used in some EOS literature to refer to only the dynamic contribution of the electron background, not the effect 
of the constant-density (static) electron background. This leads to the unfortunately worded but true statement that in highly degenerate 
plasmas where the screening contribution of the electrons to the EOS arc negligible, the screening of reaction rates by electrons can be 
quite strong. 
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and solve numerically for ji'^ and /ijf using a Newton-Raphson method. 

For charged particle reactions, rates are measured or calculated in one direction and t he reverse rate is then calculated 
based on reciprocity and the species ratio in the NSE state (see e.g. lFowler et al.lll967j) . That is, for the reaction i(j,k)l 



{(TV 



NSE 



where (ctv) denotes the product of the cross section an d velocity averaged over the appropriate distribution, fo r which 
we use tabulated values from the REACLIB database l|Thielemann TTanrmS^ lR Auscher ThieleTTiannll2nnflt) and the 
rii denote the number density of various particles. This relation continues to hold when plasma Coulomb corrections are 
important (F; ^0.1) and thus the rate enhancement factors for the forward and reverse rates must obey a relationship 
which can be derived from the NSE condition. We have foimd that traditional screening factors (e.g. I Wallace et alJ 
119821) do not inherently satisfy this relationship. Even if they did, differences in the fitting form adopted for the 
plasma Coulomb corrections would still lea d to inconsistency. Therefore we have chosen to calculate the reverse ratios 
explicitly using the same fitting form, I|A1(I . as for the NSE calculation. This leads to the relations 

7/ gkgi \mkmij \ kT J \ kT I 



((TV 



(crv)/ gi \ irik J \ J \ kT J y kT J 

for i{j,k)l and reactions respectively, where gi are the temperature dependent nuclear partition functions, and 

m,: are the nuclear masses. The forward rates themselves are found by applying screening factors from lWallace et alJ 
l)1982t) to the rates tabulated in the REACLIB database. 

Although this formalism establishes consistency between a reaction network calculation and a direct NSE calculation, 
some ambiguity remains. A favored reaction direction must be chosen and the reverse rate computed from its screened 
rate. The choice is apparent in the case of photodisintegrating reverse reactions, but in the case of {a,p) reactions, for 
example, the choice is less clear. In practice one direction typically has a superior calculation or measurement of the 
nuclear cross section, or one direction is not actually known at all. As shown below, it is not necessary to reconcile 
this ambiguity in the current work, but it is enlightening to examine the relationship between the screening corrections 
to the rates and the reverse ratio presented in (|A6(I and l|A7|l . The screening enhancement / = exp{H) is given by 
I Wallace et all ()1982() (see also lOeata et al.iri993l for a similar expression), 

H = C-^ ( J-/7^-0.014/7'^-0.128Z7^ )-F(0.0055/7'*-0.0098fo^ + 0.0048Z7^) (A8) 



3 V32 

where = 3Feff/ri2, Tn = [(27tt^ /4-)2AZjZ^e'^ /NAkTh^y^^ = 4.2(ZjzlAi2/Tgy^\ and A = AiA2/(Ai +A2) are parameters 
related to the quantum contributions and the classical plasma contributions are contained in 

C = /ex(Z^/'F,) + /ex(Z2/'F,)-/ex[(Zi +Z2)^/^F,] (A9) 

where /ex(r,) = l^f /kT is the excess free energy due to Coulomb corrections. For the situations under consideration 
here, F < 10 and Tq ^^4-9, so that b is fairly small. In this situation the classical correction, C, dominates, so that for 
an iij,k)l reaction {a\')f= (o'v)^ exp(C,j), and where B denotes the bare (unscreened) reaction rate. We thus find that 
because Z,+Zj = Zj^+Z/, the reverse ratio equation IjAGfl changes the C,j multiplying the forward rate into Q/, giving the 
the proper screening enhancement for the reverse rate. This does not address differences in the fitting form adopted 
for /ex(r(), which should be small by definition, but which we have avoided by using the same /ex(r/) in the reverse 
ratios and the NSE calculation. 
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TABLE 3 

Temperature, energy release, H and He mass fractions for isochoric, self-heating network calculations with three 

different networks. 
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No te. — P200: the 200 nuclide netw ork described i n this paper; aproxlQ: the approximated 19 nucUde network iWeaver et al.|1978i iTimme j 
119991) : torch47: the 47 nuclide network iTimmeJl99gl) . 
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TABLE 4 

Temperature, energy release, H and He mass fractions for isobaric, 
self-heating network calculations with 200 nuclide network. 
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7.80 


0932 


0.0181 


24.2 


7.10 


6.72 


3.77 


7.75 


0.149 


0.0397 


24.3 


7.18 


6.82 


3.94 


7.69 


0.230 


0826 


24.4 


7.26 


6.91 


4.13 


7.61 


0.343 


0.163 


24.5 


7.34 


7.00 


4.31 


7.49 


0.491 


0.303 


24.6 


7.41 


7.09 


4.50 


7.33 


0.674 


0.533 


24.7 


7.49 


7.18 


4.70 


7.14 


0.886 


0.882 


24.8 


7.57 


7.28 


4.89 


6.91 


1.12 


0.138 


24.9 


7.64 


7.37 


5.07 


6.66 


1.35 


2.03 


25.0 


7.72 


7.46 


5.26 


6.40 


1.57 


2.83 


25.1 


7.80 


7.55 


5.44 


6.12 


1.77 


3.77 


25.2 


7.87 


7.64 


5.61 


5.85 


1.94 


4.79 


25.3 


7.95 


7.74 


5.78 


5.60 


2.09 


5.87 


25.4 


8.02 


7.83 


5.94 


5.35 


2.21 


6.95 


25.5 


8.10 


7.91 


6.10 


5.12 


2.31 


8.01 


25.6 


8.18 


8.00 


6.25 


4.91 


2.40 


9.01 


25.7 


8.25 


8.09 


6.40 


4.72 


2.47 


9.95 


25.8 


8.33 


8.17 


6.55 


4.54 


2.53 


10.8 


25.9 


8.40 


8.26 


6.70 


4.38 


2.58 


11.6 


26.0 


8.48 


8.34 


6.84 


4.24 


2.62 


12.4 


26.1 


8.55 


8.43 


6.99 


4.11 


2.66 


13.0 


26.2 


8.63 


8.51 


7.13 


3.99 


2.69 


13.6 


26.3 


8.70 


8.59 


7.28 


3.88 


2.72 


14.2 


26.4 


8.78 


8.68 


7.43 


3.78 


2.75 


14.6 


26.5 


8.85 


8.76 


7.58 


3.69 


2.77 


15.1 


26.6 


8.93 


8.84 


7.73 


3.61 


2.80 


15.5 


26.7 


9.00 


8.92 


7.89 


3.54 


2.82 


15.8 


26.8 


9.08 


9.00 


8.04 


3.47 


2.84 


16.1 


26.9 


9.15 


9.08 


8.21 


3.41 


2.86 


16.4 


27.0 


9.23 


9.16 


8.37 


3.36 


2.88 


16.6 


27.1 


9.31 


9.24 


8.54 


3.31 


2.90 


16.8 


27.2 


9.38 


9.32 


8.72 


3.26 


2.92 


17.0 


27.3 


9.46 


9.40 


8.90 


3.23 


2.93 


17.1 


27.4 


9.53 


9.47 


9.08 


3.19 


2.95 


17.2 


27.5 


9.61 


9.55 


9.27 


3.16 


2.97 


17.2 


27.6 


9.68 


9.63 


9.47 


3.14 


2.98 


17.3 


27.7 


9.76 


9.71 


9.67 


3.11 


3.00 


17.3 


27.8 


9.83 


9.79 


9.88 


3.09 


3.02 


17.3 


27.9 


9.91 


9.86 


10.1 


3.08 


3.03 


17.2 


28.0 


9.98 


9.94 


10.3 


3.07 


3.04 


17.2 


28.1 


10.1 


10.0 


10.5 


3.06 


3.06 


17.2 


28.2 


10.1 


10.1 


10.8 


3.05 


3.06 


17.1 


28.3 


10.2 


10.2 


11.0 


3.05 


3.06 


17.0 


28.4 


10.3 


10.3 


11.3 


3.05 


3.07 


16.9 


28.5 


10.4 


10.3 


11.5 


3.05 


3.08 


16.7 



